1. Field of the Invention
The invention relates generally to imaging in a scattering medium, and more particularly to a modification of the Normalized Difference Method for solving the perturbation equation associated with a DYnamic Near-infrared Optical Tomography (DYNOT) imaging system, with improved speed and efficiency, to enable image recovery in real time.
2. Description of Related Art
Many techniques and systems have been developed to image the interior structure of a turbid medium through the measurement of energy that becomes scattered upon being introduced into a medium. Typically, a system for imaging based on scattered energy detection includes a source for directing energy into a target medium and a plurality of detectors for measuring the intensity of the scattered energy exiting the target medium at various locations with respect to the source. Based on the measured intensity of the energy exiting the target medium, it is possible to reconstruct an image representing the cross-sectional scattering and/or absorption properties of the target. Exemplary methods and systems were disclosed in Barbour et al., U.S. Pat. No. 5,137,355, entitled “Method of Imaging a Random Medium,” the disclosure of which is hereby incorporated by reference for all purposes, and Barbour, U.S. Pat. No. 6,081,322, entitled “NIR Clinical Opti-Scan System” the disclosure of which is hereby incorporated by reference for all purposes.
In particular, DYnamic Near-infrared Optical Tomography (DYNOT) has shown great promise in the field of non-invasive medical imaging. Optical tomography permits the use of near infrared energy as an imaging source. Near infrared energy is highly scattered by human tissue and is therefore an unsuitable source for transmission imaging in human tissue. However, these properties make it a superior imaging source for scattering imaging techniques. The ability to use near infrared energy as an imaging source is of particular interest in clinical medicine because the intensity of the detected energy is exceptionally responsive to blood volume and blood oxygenation levels, thus having great potential for detecting cardiovascular disease, tumors and other disease states.
More recent refinements of these techniques are disclosed in WO 01/20546, published 22 Mar. 2001, “Imaging of Scattering Media Using Relative Detector Values,” filed 14 Sep. 2000 as PCT/US00/25156, the disclosure of which is hereby incorporated by reference for all purposes. Therein, a Normalized Difference Method is described that can limit the effect of modeling errors, minimize ill-conditioning of the perturbation equation, and generally stabilize the solution to the image reconstruction problem.
In particular, it was noted that a common approach for the reconstruction of an image of the cross-sectional properties of a scattering medium is to solve a perturbation equation based on the radiation transport equation, which describes the propagation of energy through a scattering medium. The perturbation formulation relates the difference between coefficient values of the true target and a specified reference medium, weighted by a proportionality coefficient whose value depends on, among other things, the source/detector configuration and the optical properties of the medium. Tomographic measurements consider some array of measurement data, thus forming a system of linear equations having the formu−ur≡δu=Wrδx,  (1)where δu is the vector of differences between a set of measured light intensities (u) and those predicted for a selected reference medium (ur), Wr is the Jacobian operator, and δx is the position-dependent difference between one or more optical properties of the target and reference media (e.g., a change in absorption coefficient, δμa, a change in the reduced scattering coefficient δμ′s, or, in the diffusion approximation, a change in the diffusion coefficient δD, where D=1/[3(μa+μ′s)]. The operator, referred to as the weight matrix, has coefficient values that physically represent the fractional change in light intensity at the surface caused by an incremental change in the optical properties at a specified point in the medium. Mathematically, this is represented by the partial differential operator ∂ui/∂xj, where i refers to the ith source/detector pair at the surface of the medium, and j to the jth pixel or element in the medium.
Although the perturbation equation in Eq. (1) can be solved using any of a number of available inversion schemes, the accuracy and reliability of the results obtained can vary greatly due to uncertainties and errors associated with the quality of the measurement data, inaccuracies in the physical model describing light propagation in tissue, specification of an insufficiently accurate reference state, the existence of an inherently underdetermined state caused by insufficiently dense measurement sets, weak spatial gradients in the weight function, and so forth.
A matter of considerable concern is the accuracy with which the reference medium can be chosen. An accurate reference is one that closely matches the external geometry of the target medium, has the same size, nearly the same internal composition, and for which the locations of the measuring probes and their efficiency coincide well with those used in the actual measurements. While such conditions may be easily met in numerical simulations, and perhaps also in laboratory phantom studies, they represent a much greater challenge in the case of tissue studies. Confounding factors include the plasticity of tissue (which deforms upon probe contact), its mainly arbitrary external geometry and internal composition, and the considerable uncertainty stemming from the expected variable coupling efficiency of light at the tissue surface. The influence of these uncertainties can be appreciated when it is recognized that the input data vector for the standard perturbation formulation (i.e., Eq. (1)) is actually the difference between a measured and a computed quantity. This vector contains information regarding the subsurface properties of the target medium that, in principle, can be extracted provided an accurate reference medium is available.
However, there are two significant concerns that are frequently encountered in experimental studies and are not easily resolvable, especially in the case of tissue studies. One concern is the expected variable coupling efficiency of light entering and exiting tissue. Nonuniformity in the tissue surface, the presence of hair or other blemishes, its variable deformation upon contact with optical fibers, and the expected variable reactivity of the vasculature in the vicinity of the measuring probe all serve to limit the ability to accurately determine the in-coupling and out-coupling efficiencies of the penetrating energy. Variations in the coupling efficiency will be interpreted by the reconstruction methods as variations in properties of the target medium and can introduce gross distortions in the recovered images. The noted concern can be minimized by adopting absolute calibration schemes; however, the variability in tissue surface qualities will limit reliability and stability of these efforts.
A second concern stems from the underlying physics of energy transport in highly scattering media. One effect of scattering is to greatly increase the pathlength of the propagating energy. Small changes in the absorption or scattering properties of the medium can, depending on the distance separating the source and detector, greatly influence the intensity of emerging energy. This consideration has important implications for the required accuracy with which the reference medium must be specified. In the context of perturbation formulations, the reference medium serves to provide predictions of the measured energy intensities and also provides the needed weight functions that serve as the imaging operators. The difficulty is that the computed reference intensities are extremely dependent on the optical coefficient values of the reference medium. Significantly, this dependence is a nonlinear function of the distance between source and detector. It follows that a small change in the optical properties of the reference medium may influence the value of the computed intensity differences (δu) by a relative amount that may be significantly different for each source/detector pair, thereby altering the information content of the data vectors. This can lead to the recovery of grossly corrupted images. While in principle such effects may be overcome by use of recursive solutions to the perturbation equation (i.e., Newton-type updates), in practice this can require extensive computational efforts, especially in the case of 3D solutions. Moreover, such efforts to improve on first-order solutions to the perturbation equation (e.g., Born or Rytov solutions), can fail if the initial estimate chosen for the reference medium is insufficiently accurate.
WO 01/20546 addressed the need for image reconstruction techniques based on the detection of scattered energy that (1) do not require absolute calibration of, and absolute measurements by, the detectors and other elements of the apparatus, (2) make the standard perturbation equation less susceptible to variations between boundary conditions and properties of the reference medium and the target medium, and (3) produce solutions having physical units.
It was noted that the typical process for imaging of a scattering medium comprises: (1) selecting a reference medium having known boundary conditions and optical properties which are substantially similar to those of the intended target; (2) determining a weight matrix and an intensity of emerging energy exiting the reference medium at each of a plurality of source points for each of a plurality of detectors located around the reference medium boundary, the determination being made by either actual measurements or solution of the radiation transport equation or an approximation thereto; (3) measuring actual emerging energy intensities for corresponding source and detector points on a target medium; and (4) solving the perturbation equation for the optical properties of the target based on the measured intensities of energy emerging from the target.
WO 01/20546 described an improved methodology for imaging of a scattering medium using a modified form of the standard perturbation equation, which was capable of (1) reducing the sensitivity of the perturbation equation to differences between the reference medium and target medium, (2) producing solutions to the perturbation equation having physical units, and (3) reducing the effect of variable detector efficiencies without the need for absolute calibration, while at the same time preserving the ability to compute recursive solutions. Compared to the standard perturbation approach, the described approach provided a remarkable improvement in the quality of image reconstruction.
While the method was applicable to known static imaging techniques, it was also instrumental in the realization of practical dynamic imaging of highly scattering media. A first element of dynamic imaging is the development of a fast, parallel, multi-channel acquisition system that employs geometrically adaptive measuring heads. This system is described in further detail in PCT/US00/25155, published on Mar. 22, 2001 as WO 01/20306, incorporated herein by reference. The second element is to evaluate the acquired tomographic data using the modified perturbation methods of the present invention. The third element is to collect a time series of data and subject either the time series of data or a time series of images reconstructed from the data to analysis using various linear and nonlinear time-series analysis methods to extract dynamic information and isolated dynamic information. These methods are described in further detail in PCT/US00/25136, published on Mar. 22, 2001 as WO 01/20305, incorporated herein by reference.
Some of the methods, systems and experimental results described focused on optical tomography of human tissue using wavelengths in the near infrared region for the imaging source. The techniques were generally applicable to the use of essentially any energy source (e.g., electromagnetic, acoustic, and the like) on any scattering medium (e.g., body tissues, oceans, foggy atmospheres, earth strata, industrial materials) so long as diffusive-type mechanisms are the principal means for energy transport through the medium.
Numerous imaging systems such as those disclosed in the previously mentioned Barbour '355 patent and the Barbour '322 patent have been developed for use in imaging of a scattering medium. A schematic illustration of an exemplary system is shown in FIG. 1. This system includes a computer 102, sources 104, 106, a source demultiplexer 108, an imaging head 110, detectors 112 and a data acquisition board 114. A target 116 placed in the imaging head 110 is exposed to optical energy from sources 104,106. The optical energy originating from sources 104,106, is combined by beam splitter 118 and is delivered to source demultiplexer 108. The source demultiplexer 108 is controlled by computer 102 to direct the optical energy to source fibers 120 sequentially.
Each source fiber 120 carries the optical energy from the demultiplexer 108 to the imaging head 110, where the optical energy is directed into the target 116. The imaging head 110 contains a plurality of source fibers 120 and detector fibers 122 for transmitting and receiving light energy, respectively. Each source fiber 120 forms a source-detector pair with each detector fiber 122 in the imaging head 110 to create a plurality of source-detector pairs. The optical energy entering the target 116 at one location is scattered and may emerge at any location around the target 116. The emerging optical energy is collected by detector fibers 122 mounted in the imaging head 110.
The detector fibers 122 carry the emerging energy to detectors 112, such as photodiodes or a CCD array, that measure the intensity of the optical energy and deliver a corresponding signal to a data acquisition board 114. The data acquisition board 114, in turn, delivers the data to computer 102. The computer 102 processes values of the measured energy for use in solving the system equation for the target medium. This imaging process is repeated so as to deliver optical energy to each of the source fibers sequentially, a measurement being obtained for detected emerging energy at each detector for each emitting source fiber. This process may continue over a period of time with the computer 102 storing the data for reconstruction of one or more images. Additionally, the system may include two or more imaging heads for comparing one target to another. The computer 102 reconstructs an image representative of the internal optical properties of the target by solving a perturbation equation. It will be appreciated by those skilled in the art that more than one computer can be used to increase data handling and image processing speeds.
The Standard Perturbation Formulation
As discussed above, reconstruction of a cross-sectional image of the absorption and/or scattering properties of the target medium is based on the solution of a perturbation formulation of either the radiation transport equation or the diffusion equation. The perturbation method assumes that the composition of the unknown target medium deviates by only a small amount from a known reference medium. This reduces a highly nonlinear problem to one that is linear with respect to the difference in absorption and scattering properties between the target medium under investigation and the reference medium. The resulting standard perturbation equation has the form of Eq. (1), discussed above, where δu is a vector of source-detector pair intensity differences between the measured target medium and the known reference medium (i.e., δu=u−ur). W is the weight matrix describing the influence that each volume element (“voxel”) of the reference medium has on energy traveling from each source to each detector, for all source-detector pairs. The volume elements are formed by dividing a slice of the reference medium into an imaginary grid of contiguous, non-overlapping pieces or sections. Physically, the weight matrix contains the first-order partial derivatives of the detector responses with respect to the optical coefficients of each volume element of the reference medium. δx is the vector of differences between the known optical properties (e.g., absorption and scattering coefficients) of each volume element of the reference medium and the corresponding unknown optical properties of each volume element of the target medium.
This standard perturbation equation assumes (1) use of absolute detector measurements for u, and (2) that the boundary conditions and optical properties of the reference do not vary significantly from those of the target. Both of these factors are problematic in practice.
The Modified Perturbation Formulation
WO 01/20546 discussed modifying the standard perturbation equation by replacing δu with a proportionate relative difference between two measured values multiplied by a reference term with the required units as set forth in the equation (2) below:
                                          (                          δ              ⁢                                                          ⁢                              I                r                                      )                    i                =                              [                                                            I                  i                                -                                                      (                                          I                      0                                        )                                    i                                                                              (                                      I                    0                                    )                                i                                      ]                    ⁢                                    (                              I                r                            )                        i                                              (        2        )            where i is the source/detector pair index. In equation (2), Ir is the computed detector reading corresponding to a source-detector pair of a selected reference medium, and I and I0 represent two data measurements for a corresponding source-detector pair on one or more targets (e.g., background vs. target, or time-averaged mean vs. a specific time point, etc.). The resultant term δIr therefore represents a relative difference between two sets of measured data that is then mapped onto a reference medium. This modification has important attributes that limit the effects of modeling errors and minimize ill-conditioning of the inverse problem (i.e., image reconstruction problem), while retaining the correct units in the solution.
The corresponding perturbation equation using this modified term is:Wr·δx=δIr  (3)In equation (3) Wr and δx are the same as Wr and δu in equation (1). Referring to equations (2) and (3), it can be seen that in the limit, when Ir equals to I0, this equation reduces to the standard perturbation formulation shown in equation (1). This formulation represents the Born approximation formulation of the modified perturbation equation. A similar expression may be written for the Rytov approximation in the following form:
                                                                                                              (                                          δ                      ⁢                                                                                          ⁢                                              I                        ′                                                              )                                    i                                =                                  ln                  ⁢                                                            I                      i                                                              (                                              I                        0                                            )                                                                                  ;                                                                                                            (                                      W                    r                    ′                                    )                                =                                                                            (                                              W                        r                                            )                                        ij                                                                              (                                              I                        r                                            )                                        i                                                              ;                                                                                          δ                ⁢                                                                  ⁢                                  I                  ′                                            =                                                W                  r                  ′                                ⁢                                                                  ⁢                δ                ⁢                                                                  ⁢                x                                                                        (        4        )            
Equation (2) preserves the information content of a measured proportionate relative data vector obtained from the target medium and yields data vector having the correct physical units. Apart from simplifying measurement requirements, the method represented by equations (3) and (4) also reduces susceptibility of the perturbation equation to boundary and optical property variations between the target and the reference media. Additionally, iterative solutions of equations (3) and (4) can be easily implemented. In principle, the techniques used in the modified perturbation equation, referred to as the normalized difference method (NDM), can be used to evaluate any measured proportionate relative quantity.
Experimental Validation
WO 01/20546 presented examples that illustrate the benefits of applying the NDM approach for the analysis of relative measures from highly scattering media.
Forward Model and Data Acquisition Geometry
For any intended target, the perturbation equation is generated for a reference medium having boundary conditions and optical properties substantially similar to the target. The perturbation equation models the energy propagation, e.g., light, in the reference medium as a diffusion process. For a domain Ω having a boundary δΩ this is represented by the expression∇·[D(r)∇u(r)]−μa(r)u(r)=δ(r−rs), rεΩ  (5)where u(r) is the photon density at position r, rs is the position of a DC point source, and D(r) and μa(r) are the position-dependent diffusion and absorption coefficients, respectively. The diffusion coefficient is defined asD=1/{3[μa(r)+μ′s(r)]}  (6)where μ′s(r) is the reduced scattering coefficient. The photon density values at the detectors, i.e., the calculated intensity of energy emerging from the reference medium at each detector, were computed by applying Dirichlet boundary conditions on an extrapolated boundary. Depending on the target medium to be explored, sources and detectors for the reference are positioned 1 to 2 transport mean free pathlengths within the boundary of the reference medium.
Solutions to the diffusion equation may be computed by any known means, such as by the KASKADE adaptive finite element method [R. Beck, R. Erdmann and R. Roitzsch, “Kaskade 3.0—An object-oriented adaptive finite element code,” Technical report TR 95-4, Konrad-Zuse-Zentrum fur Informationstechnik, Berlin (1995)]. This is a publicly available code suitable for the solution of partial differential equations in one, two or three dimensions using adaptive finite element techniques. The code can be readily modified to permit solutions to the diffusion equation using a point source. Mesh generation may be by any known method, such as the Delaunay tessellation algorithm originally described by Watson [D. F. Watson, “Computing the n-dimensional Delaunay tessellation with applications to Voronoi polytopes”, Computer Journal, 24,167-172 (1981)].
The perturbation equation is specific to the boundary conditions and optical properties of the reference medium, including the orientation of the source-detector pairs in relation to one another and to the reference medium. These conditions and properties are preferably nearly identical to those of the target. For example, in the experiments discussed below, the perturbation equation was generated based on an imaging system having six sources and eighteen detectors per source (108 source-detector pairs) with the sources equally spaced at 60 degree intervals around the boundary of the medium and the detectors equally spaced at 20 degree intervals.
Inverse Algorithm
As described above, the relative intensities are measured for all source-detector pairs using any known imaging system. Image recovery is then achieved using known methods, such as conjugate gradient descent (CGD), or simultaneous algebraic reconstruction techniques (SART), to solve the modified perturbation equation for the absorption and scattering properties of the target. See J. Chang, H. L. Graber, R. L. Barbour and R. Aronson, “Recovery of optical cross-section perturbations in dense-scattering media by transport-theory-based imaging operators and steady-state simulated data”, Appl. Opt. 35, 3963-3978, (1996) (the disclosure of which is incorporated herein by reference). For example, the experimental results discussed below were achieved using a CGD solver, both with and without matrix rescaling. In addition, a weight matrix resealing (WMR) technique may be used to improve the conditioning of the weight matrix. The effect of rescaling the weight matrix is to make it more uniform. Two rescaling criteria can be applied for this purpose: (1) rescaling the maximum of each column to 1; or (2) rescaling the average of each column to 1. When WMR was used, criterion (1) was applied for image recovery. The solution to the modified perturbation equation provides a relative measure of the difference between the cross-sectional optical properties of a target during the first and second measurements I and I0. The values from this solution are used to generate cross-sectional images representative of the target's internal optical properties.
Building upon the above knowledge, the present invention modifies the Normalized Difference Method to enable real-time image recovery in a tomography or other imaging system. Because of the inherent stability of solutions obtained from the Normalized Difference Method, a variety of computation-saving techniques and stability-enhancing techniques are introduced to provide fast, accurate reconstruction of a time series of images involving large-scale 3D problems.